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(57) ABSTRACT 

A system for simulating operation of a VLSI interconnect 
structure having capacitive and inductive coupling between 
nodes thereof, including a processor, and a memory, the pro- 
cessor configured to perform obtaining a matrix X and a 
matrix Y containing different combinations of passive circuit 
element values for the interconnect structure, the element 
values for each matrix including inductance L and inverse 
capacitance P, obtaining an adjacency matrix A associated 
with the interconnect structure, storing the matrices X, Y, and 
A in the memory, and performing numerical integration to 
solve first and second equations. 
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COMPUTATIONALLY EFFICIENT 
MODELING AND SIMULATION OF LARGE 
SCALE SYSTEMS 

CROSS-REFERENCE TO RELATED 
APPLICATIONS 


Gx 4 - Cx = b. 


where 


Q 

-Ai 




(1) 


This application is a continuation application of patent 
application Ser. No. 12/852,942, filed on Aug. 9, 2010, now 
U.S. Pat. No. 8,336,014 to Jain et al., issued Dec. 18, 2012, 
which is a divisional application of patent application Ser. No . 
1 1/593,465, filed on Nov. 6, 2006, now U.S. Pat. No. 7,774, 
725 to Jain et al., issued Aug. 10, 2010, which claims the 
benefit of Provisional Patent Application No. 60/733,460, 
filed Nov. 4, 2005, and Provisional Patent Application No. 
60/740,990, filed Nov. 30, 2005, which applications are 
hereby incorporated by reference along with all references 
cited therein. 


GOVERNMENT RIGHTS 

This invention was made with government support under 
Contract/Grant No. NCC 2-1363 awarded by the National 
Aeronautics and Space Administration (NASA), under Con- 
tract/Grant Nos. CCR-9984553 and CCR-0203362 awarded 
by the National Science Foundation, and under Contract/ 
Grant No. USAF-FA8650-04-D-2409 awarded by the United 
States Air Force Research Laboratories. The government has 
certain rights in the invention. 
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, Q = A T R~ l A g , and C = A T c CA c . 
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R is the resistance matrix. The matrices Q , L and C are the 
conductance, inductance and capacitance matrices respec- 
tively, with corresponding adjacency matrices A^, A I and A c . 
15 \ s i s the current source vector with adj acency matrix A,., and v n 
and i ; are the node voltages and inductor currents respectively. 
With n denoting the number of inductors, we note that 
L, C, Q eR 2nx2n . 

A standard algorithm for the numerical integration of dif- 
20 ferential equations such as (1) is the trapezoidal method. 
Consider a uniform discretization of the time axis with reso- 
lution h. Then, using the notation x*=x(kh), and the approxi- 
mations 


25 



x k+l -x k 

h 


and x? 1 as 


** +1 -** 

2 


30 over the interval [kh,(k+l )h] , we may solve for x* +1 in terms 
of x k by solving the equation 


TECHNICAL FIELD OF THE INVENTION 


The present invention relates generally to electrical circuit 
modeling and simulation techniques and, more particularly, 35 
to methods for simulating interconnect effects in very large 
scale integrated circuits. 

BACKGROUND OF THE INVENTION 

40 

With aggressive technology scaling, the accurate and effi- 
cient modeling and simulation of interconnect effects has 
become (and continues to be) a problem of central impor- 
tance. In a three-dimensional interconnect structure there can 
be significant amounts of coupling, both inductive and 45 
capacitive, between interconnects. Models that capture these 
effects tend to involve large matrices, resulting in extraordi- 
nary demands on memory. Simulation with these models 
require prohibitive amounts of computation. 

While all coupling effects in theory extend without bound, 50 
it is well-recognized that, in practice, the effects of capacitive 
coupling, and to some extent that of inductive coupling, can 
be assumed to be local without much sacrifice in accuracy. 
Practical modeling and simulation techniques exploit this 
localization to significantly reduce storage and computa- 55 
tional costs. For practical interconnect structures, the capaci- 
tance matrix C and the inverse of the inductance matrix 
K=L -1 turn out to be (approximately) sparse. A number of 
techniques exploit the sparsity in K at extraction level. 
Exploiting sparsity of C and K in simulation however, i s much 60 
less straightforward. The main problem is that simulation 
requires terms that not only involve the sparsified matrices C 
and K, but also inverses of terms that involve them; these 
inverses are dense in general. 

The Modified Nodal Analysis (MNA) of interconnect 65 
structures such as the one shown in FIG. 1 yields equations of 
the form 



c) . M +1 +^ 
7Tf + ^— 


(2) 


A direct implementation of this algorithm requires 0(n 3 +pn 2 ) 
operations, where p is the number of time steps. The direct 
implementation ignores the structure of the matrices G and C 
that is evident in (1); explicitly recognizing this structure 
yields the so-called Nodal Analysis (NA) equations, used in 
INDUCTWISE: 


( @+ \ c+ fY" 1 = Y e + \ c - YY ■ 2aJ + Ar</ ‘ +1 + (3) 

v 

and 

2 A] tl +i =2 Af$ + hS(v * +1 + vj), (4) 


where S=A / KA / r (recall that K=L -1 , L being the inductance 
matrix, with corresponding adjacency matrix A z , and Af 
being the transpose of A z ). 

The NA equations (3) and (4) enjoy several advantages 
over the MNA equations (1). The first advantage is that the 
solution of equations (1), a problem of size 3n, has been 
divided into two sub-problems of sizes 2n and 2n, which 
yields computational savings with polynomial-time algo- 
rithms. Next, it has been observed that with typical VLSI 
interconnect structures, the matrices K, C and Q exhibit spar- 
sity. This can be used at the extraction stage to write down (3) 
and (4) with fewer parameters. Finally, at the simulation 
stage, the structure of the matrix U defined in (3) — symmetry, 
positive-definiteness and sparsity — lends itself to the use of 
fast and sound numerical techniques such as Cholesky fac- 
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torizations. These advantages have been extensively used to 
develop INDUCTWISE. For future reference, we note that 
the computation with INDUCTWISE is 0(n 3 +pn 2 ) opera- 
tions, and is usually dominated by 0(pn 2 ). 

5 

SUMMARY OF THE INVENTION 

The approach that is employed is to sparsity the various 
matrices that underlie the model of interconnects; the result- 
ing approximate models can be represented by far fewer 10 
parameters, leading to savings in storage. 

The present invention presents methods that systematically 
take advantage of sparsity in C and K, in simulation, achiev- 
ing a very significant reduction in computation with very little 
sacrifice in simulation accuracy. The first idea underlying our 
approach is that if the sparsity in the inverse of a dense matrix 
is known, the (sparse) inverse can be computed very effi- 
ciently. We take advantage of this fact by writing the simula- 
tion equations in terms of L and P=C -1 . The most computa- 2 o 
tionally intensive step in simulation, of system formulated in 
such a fashion, reduces to that of matrix-vector multiplication 
involving a sparse matrix. We also show that savings with 
sparse-matrix-vector multiplication can be obtained with 
simulation using K=L -1 and C, as well, but to a lesser extent. 25 

The RLP formulation is extended to include non-linear 
devices, without sacrificing the computational benefits 
achieved due to sparsity of the linear system. It should be 
noted that the A matrix involved in the solution of the linear 
system is constant throughout the simulation. In contrast, the 30 
A matrix involved in solving the non-linear system changes in 
each simulation step. However, the A matrix is sparse. Due to 
the sparse and time varying nature of the problem at hand 
Krylov subspace based iterative methods could be used for 
efficient simulation. Our second contribution is to introduce a 
novel preconditioner constructed based on the sparsity struc- 
ture of the non-linear system. The inverse of the precondi- 
tioner has a compact representation in the form of the Had- 
amard product, which facilitates not only the fast 40 
computation of the inverse, but also the fast dense matrix- 
vector product. 

The objects and advantages of the present invention will be 
more apparent upon reading the following detailed descrip- 
tion in conjunction with the accompanying drawings. 45 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 a distributed model of a typical three dimensional 
VLSI interconnect structure that may be operated in Simula- 50 
tion with the methods of the present invention. 

FIG. 2 shows the average sparsity index of the matrices 
U -1 V, U _1 A z r , U -1 A z r and U -1 S, for a structure as a fimction 
of h for various values of €. 

FIG. 3 shows average sparsity index of the matrices X -1 Y, 55 
X -1 A and X _1 AP, for a structure as a fimction of h for the 
sparsity threshold of €=0.001, as compared with the average 
sparsity index of the matrices encountered in the GKC-algo- 
rithm. 

FIG. 4(a) shows the significant entries (shown darker) of 60 
the adjacency matrix A for a structure with 1500 conductors. 

FIGS. 4(b) and 4(c) show the significant entries (shown 
darker) of W _1 andX -1 , respectively, for the structure in FIG. 
4(a ) . 

FIG. 5 shows the voltage wave forms, obtained from 65 
SPICE and Exact-RLP, of the active line and the seventh line 
of a 100-conductor circuit. 


4 

FIG. 6 shows the voltage wave forms, obtained through 
INDUCTWISE, Exact-RLP, RLP, andGKC, of the active line 
and the seventh line of a 600-conductor circuit. 

FIGS. 7 and 8 show plots of the RMSE for the active and 
the seventh line as a function of threshold value for a 600- 
conductor circuit. 

FIG. 9 shows the sparsity structure (nonzero entries shown 
darker) of the A matrix for an exemplary circuit of parallel 
wires driving a bank of inverters. 

FIG. 10 shows an exemplary preconditioner matrix that 
may be used with the exemplary circuit of FIG. 9. 

FIG. 11 shows the sparsity pattern (nonzero entries shown 
darker) of matrix A of a circuit having only non-linear devices 
and no interconnects. 

FIG. 12 shows average sparsity versus circuit size. 

FIG. 13 shows the voltage wave form obtained through 
SPICE and Exact-RLP and SASIMI. 

FIG. 14 shows a two dimensional non-uniform spatial grid 
for a Nanotransistor. 

FIG. 15 shows the ratio of memory consumption of the 
algorithm in [9] as compared to ours for varying division sizes 

oyD). 

FIG. 16 shows the ratio of memory consumption of the 
algorithm in [9] as compared to ours for a varying number of 
divisions (D). 

DETAILED DESCRIPTION OF THE PREFERRED 
EMBODIMENTS 

For the purpose of promoting an understanding of the 
principles of the invention, reference will now be made to the 
embodiments illustrated in the drawings and specific lan- 
guage will be used to describe the same. It will nevertheless 
be understood that no limitation of the scope of the invention 
is thereby intended, such alterations and further modifica- 
tions in the illustrated device and such further applications of 
the principles of the invention as illustrated therein being 
contemplated as would normally occur to one skilled in the art 
to which the invention relates. 

While significant storage and computational advantages 
accrue with INDUCTWISE, we note that the sparsity of U has 
not been fully taken advantage of at the level of linear algebra 
(beyond the possible use of sparse Cholesky factorizations) in 
the numerical solution of (3). In particular, with the formula- 
tion used by INDUCTWISE, while the matrix U is sparse, its 
inverse is dense. Thus, trapezoidal numerical integration, at a 
first glance, entails matrix-vector multiplies with a dense 
matrix at each t im e step. However, it has been observed that 
the matrices U -1 V (where V is defined in (3)), U _1 A / r , 
U -1 A z r and U -1 S are approximately sparse, and this informa- 
tion can be used to significantly reduce the computation as 
follows. Rewrite (3) and (4) as 

v„ k+ l =lT 1 Vv n k -2 U~ l A 1 r i j k + IT l A i r (I s k+l +lJ c ) 

2LT l A fi l k+i =2 LT l A ji f+h IT l S(v„ k+i +v„ k ) . 

Pre-compute and store the sparsified matrices U -1 V, 
U -1 A z r , U -1 A z r and U -1 S. Then, every time step in the trap- 
ezoidal integration scheme requires only sparse matrix-vec- 
tor multiplies. We will henceforth refer to this technique as 
the GKC-algorithm (as the computations are done with the 
conductance, inverse of the inductance and the capacitance as 
the parameters). 

In order to quantity the computational savings obtained 
with the GKC-algorithm over INDUCTWISE, we define the 
“ sparsity index” p e (A) of a matrix A as ratio of the number of 
entries of A with absolute value less than € to the total number 
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of entries. Then, the computation required for each iteration 
with the GKC-algorithm, with some appropriate value of €, is 
0((1 -v)n 2 ) where v is the minimum of the sparsity indices of 
the matrices U _1 V, U~ 1 A l T , U _1 A z r andU _1 S. The value of v 
can be expected to depend on the threshold for detecting 5 
sparsity €, as well as the time step size h. FIG. 2 shows the 
average sparsity index of the matrices U _1 V, U -1 Af, U -1 A z r 
and U _1 S, for a structure with three parallel planes consisting 
of 600 conductors, as a function of h for various values of €. 
The typical value of h used solving the MNA equations for 10 
VLSI interconnects is 0. 1 picoseconds. With such values of h 
and €=0.001 , it can be seen that v»0.8. Thus the total compu- 
tation time with the GKC-algorithm is approximately a fifth 
of that required by INDUCTWISE. 15 

We now explore an alternative formulation of the MNA 
equations that uses the resistance, inductance and the inverse 
of the capacitance matrix. For typical interconnect structures, 
shown in FIG. 1, we can manipulate the MNA equations (2) to 
obtain 


( L R h . 

U + 2 + 4 APA T = 


(5) 


(L 

R 

h T \ , 

- - 



- -apa t ; 

u 

2 

4 1 


and 


„*+! - v * _ -PA r $ + 1 + if) + -P(I* +l + /*), 


( 6 ) 30 


where P is the inverse capacitance matrix i.e P=C -1 and A is 
the adj acency matrix of the circuit, obtained by first adding A^ 3 
and A 1 and then removing zero columns (these correspond to 
intermediate nodes, representing the connection of a resis- 
tance to an inductance). When compared with the NA equa- 
tions (3) and (4), we see that the number of state variables has 
been halved. Compared to INDUCTWISE, this represents 4Q 
immediate savings. For future reference, we will term the 
technique of directly solving (5) and (6) as the “Exact-RLP” 
algorithm. 

In contrast with the GKC-algorithm, it turns out here that X 
is dense, but with an inverse that is approximately sparse. 
Thus, windowing techniques such as those employed by 45 
INDUCTWISE during the extraction stage to obtain a spar- 
sified matrix K can be employed here to quickly compute a 
sparsified X -1 . (Windowing techniques details will be 
described below.) Moreover, the matrices X -1 Y, X -1 A and 
X -1 AP turn out to be approximately sparse. Thus, paralleling 50 
the development of the GKC-algorithm, we have the follow- 
ing RLP-algorithm: 

Rewrite (5) and (6) as 


/f +1 = x~ 1 y$ +x~ l Av k n + ^x~ l AP(i* +l + /*), 

X _ 1 Avf +1 = X _ 1 AvJ- ^X _ 1 APA r (i f +1 + if ) + ^X~ l AP(^ +l + /*). 

60 

Pre-compute and store the sparsified matrices X -1 Y, X -1 A 
and X -1 AP. Again, every time-step in the trapezoidal integra- 
tion scheme requires only sparse matrix-vector multiplies. As 
with the GKC-algorithm, the total computation with the RLP- 
algorithm is dominated by 0((l-y)n 2 ), where is the y is the 65 
minimum of the sparsity indices the matrices X -1 Y, X -1 A and 
X _1 AP. 


FIG. 3 shows the average sparsity index of the matrices 
X -1 Y, X -1 A and X _1 AP, for a structure with three parallel 
planes consisting of 600 conductors, as a function of h for the 
sparsity threshold of €=0.001, as compared with the average 
sparsity index of the matrices encountered in the GKC-algo- 
rithm. It is clear that the matrices encountered in the RLP- 
algorithm exhibit much higher sparsity over a wide range of 
time- steps. In particular, for h=0.1 ps, it can be seen that 
y«0.9. Thus the total computation time with the above RLP- 
algorithm is approximately one-tenth of that required by the 
RLP formulation that does not use sparsity information. 
When compared to the GKC-algorithm and INDUCTWISE 
which use twice as many state variables, the amount of com- 
putation required by the RLP-algorithm is approximately 
one-eighth and one-fiftieth respectively. 

We now provide the details on the fast inversion of X. 
Assume for simplicity that the sparsity pattern in X -1 is 
known, deferring for later the problem of detecting this spar- 
sity pattern. Then, manipulations of only a subset of the 
entries of the X matrix (rather than the entire X matrix) can be 
used to compute the inverse matrix. To briefly illustrate the 
idea consider the example when XeR 1 1 x 1 1 , and the 5th row of 
X -1 has the following form: 

[0 0 ★ 0 ★ 0 0 ★ 0 0 ], 

where ★ denotes the nonzero entries. Then, it can be shown 
that these nonzero entries can be computed exactly from the 
second row of the inverse of the following 3x3 matrix 
obtained from X: 


x 33 x 35 x 38 - 

X 5 3 X 55 X 58 . 

x 83 X 8 s x 88 


More generally, suppose that there are a z nonzero entries in 
the ith row of X -1 . By following a procedure as above, the ith 
row of X -1 can be computed by inverting an a z xa z matrix. 
Thus, the overall computation in determining is X -1 is 
0(2 z a z 3 ). It is typical with VLSI interconnects that a z is a 
small constant. Thus if X -1 is exactly sparse, with a known 
sparsity pattern, it can be computed in O(n) from X. Table 1 
gives the time taken for inversion for different circuit sizes. 


TABLE 1 


Inversion time in matlab (in seconds) 



No. 

of conductors 



500 

1000 

2000 

5000 

Direct Inversion 

.29 

2.18 

16.87 

260.68 

Fast Inversion 

.79 

1.48 

2.93 

10.17 


Thus, there remains the problem of determining the spar- 
sity pattern in X -1 . Recall that 


L R 


-apa t . 


Let 


L R h 

- + — and Z = -A PA 
hi 4 


T 
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Then 

x- 1 =fr 1 -fr 1 (fr 1 +z- 1 )- 1 fr 1 . (7) 

For the values of R, L, C and h under consideration, it turns 
out that 

xr^vr'-vr'zjrK ( 8 ) 

Thus, the significant entries of X -1 can be obtained by super- 
posing the significant entries of W -1 and the significant 
entries of W -1 ZW -1 . The sparsity pattern of W _1 can be 
efficiently determined using the techniques available in the 
literature. Turning next to 

W~ l ZW~ l = - A W~ l APA T W ~ l , 

4 

note that the significant entries of W -1 A are obtained by 
distributing the significant entries of W -1 into locations deter- 
mined by the adjacency matrix A. In summary, we have the 
following heuristic for predicting the sparsity pattern in X -1 : 
First determine the significant entries of W _1 by determining 
the set of segments that are inductively couple with a given 
segment. In addition, spread the nonzero entries of W -1 to 
locations suggested by the adjacency matrix to find the 
remaining significant entries. 

These ideas are illustrated via a three dimensional inter- 
connect structure of three parallel planes with 1500 conduc- 
tors. In FIG. 4 (a), the significant entries of the adjacency 
matrix A are shown to be darker. FIGS. 4(b) and 4(c) show the 
entries ofW -1 andX -1 respectively, again with the significant 
entries shown darker. 

We emphasize that the actual computation of the signifi- 
cant entries of X -1 proceeds via the technique in, where given 
the knowledge of the sparsity pattern resident in X” 1 , the 
actual entries can be directly and efficiently computed. Thus, 
(7) and (8) are not used for computation, but only to motivate 
the heuristic for efficiently determining the sparsity pattern of 
X" 1 . 

We implemented the INDUCTWISE, RLP and GKC algo- 
rithms in MATLAB on a PC with an Intel Pentium IV 2 .4 GHz 
processor. In order to quantity the simulation accuracy with 
various methods, we used as the benchmark the Exact-RLP 
simulation (recall that this is the direct simulation of equa- 
tions (5) and (6)). (While SPICE simulations would have been 
more natural to use as the benchmark, we found that the 
computation time grew quickly to make them impractical; for 
a modest-size circuit comprising 100 parallel conductors, 
SPICE simulation took 350 seconds as compared to 1.08 
seconds with the Exact-RLP algorithm, with no detectable 
simulation error, as shown in the FIG. 5). 

Simulations were done on a three dimensional structure of 
three parallel planes, with each plane consisting of busses 
with parallel conductors, with wire-lengths of 1 mm, and a 
cross section of 1 jamx 1 pm. The wire separation was taken to 
be 1 pm; each wire was divided into ten segments. A periodic 
IV square wave with rise and fall times of 6 ps each was 


8 

applied to the first signal on the lowest plane, with a time 
period of 240 ps. All the other lines were assumed to be quiet. 
For each wire, the drive resistance was 10Q the load capaci- 
tance was 20 fF. A time step of 0.15 ps was taken and the 
5 simulation was performed over 330 ps (or 2200 time steps). 
As expected, with all methods, there is an inherent trade- 
off between simulation accuracy and cost (CPU time and 
memory). We first present results comparing the accuracy in 
simulating the voltage waveforms at the far end of the first 
to (active) and the seventh (victim or quiet) lines. The metric for 
comparing the simulations is the relative mean square error 
(RMSE) defined as 

15 J 

Z v? 

where v and v denote the waveforms obtained from Exact- 
RLP and the algorithm under consideration respectively. A 
20 threshold value of 0.001 was chosen for sparsification of RLP 
and GKC algorithms, as well as for sparsification of L -1 in 
INDUCTWISE. Table 2 presents a summary of the results 
from the study of simulation accuracy. 

25 TABLE 2 




RMSE 

comparison. 





Active Line 



7th line 


Size 

INDUCTWISE 

RLP 

GKC INDUCTWISE RLP 

GKC 

300 

.0013 

.0010 

.0017 

.1622 

.1266 

.1960 

600 

.0014 

.0011 

.0014 

.4381 

.3452 

.4651 

900 

.0006 

.0007 

.0008 

.3222 

.3076 

.4078 

1200 

.0004 

.0004 

.0004 

.2382 

.2656 

.2992 

1500 

.0003 

.0003 

.0004 

.2021 

.2200 

.2336 


It can be seen that the simulation accuracy of the RLP and 
the GKC algorithms are comparable to that of INDUCT- 
WISE, with a marginally inferior performance as measured 
40 by the RMSE. A plot of the voltage waveforms at the far end 
of the active line and the 7th line, obtained by INDUCT- 
WISE, RLP, and GKC algorithms, is shown in the FIG. 6. 

We briefly explore the influence the choice of the threshold 
for determining sparsity. A higher threshold can be expected 
45 to decrease the computational and memory requirements, 
however with loss in simulation accuracy. FIGS. 7 and 8 show 
plots of the RMSE for the active and seventh line as a function 
of threshold value, again for a circuit of size 600 conductors. 
Any value of the threshold below 0.001 appears to be a 
50 reasonable choice. 

We now turn to a comparison of the computational and 
memory requirements between INDUCTWISE, RLP and 
GKC algorithms. Table 3 summarizes the findings. 


TABLE 3 




Run time and 

memory comparisons 






Time (in sec) 




Memory (in MB) 



Size 

Exact-RLP 

INDUCTWISE 

RLP 

GKC 

Exact-RLP 

INDUCTWISE 

RLP 

GKC 

300 

14.30 

74.34 

4.09 

18.99 

2.95 

11.61 

1.02 

6.61 

600 

76.21 

422.00 

16.28 

77.32 

11.61 

46.20 

2.36 

15.38 

900 

244.14 

1133.40 

33.21 

162.08 

26.03 

103.84 

4.09 

31.68 

1200 

513.08 

3051.10 

60.53 

312.93 

46.20 

184.56 

6.16 

52.22 

1500 

827.50 

4682.00 

92.16 

813.00 

72.14 

288.24 

7.60 

86.00 
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It can be seen that for a circuit consisting of 1200 conduc- 
tors, RLP is about nine times faster than the Exact-RLP, and 
fifty times faster than INDUCTWISE. The GKC algorithm is 
about twice as fast as the Exact-RLP, and ten times faster than 
INDUCTWISE. The Exact-RLP is about six times as fast as 5 
INDUCTWISE. With larger circuit sizes, the advantage of 
RLP over INDUCTWISE continues to grow, while the Exact- 
RLP and GKC algorithms have an advantage over INDUCT- 
WISE that grows slightly. An explanation for the slower per- 
formance of INDUCTWISE compared to Exact-RLP is that 10 
the number of variables with the latter algorithm is one-half as 
that with the former. The same trends are observed with 
memory requirements. 

VLSI interconnect structures with non linear devices can 
also be analyzed using the Modified Nodal Analysis (MNA) 15 
formulation, yielding equations of the form 


Gx + Cx = b, 
where 

Q A] 


G = 


-A 


0 


, C 0] 

,C= ,* = 

1 1 0 L U 


(9) 


20 


b = 


Ajls + Inl 
0 


, Q = AlR-'Ag, and C = A T c CA c . 


R denotes the resistance matrix. The matrices Q , L and C are 
the conductance, inductance and capacitance matrices 30 
respectively, with corresponding adjacency matrices A^, A l 
and A c . l s is the current source vector with adjacency matrix 
A z , and v n and i I are the node voltages and inductor currents 
respectively. 

Vector, \ nl captures the effect of non-linear loads and 35 
depends on the node voltages as I w/ =f(v w ). f is a function 
which varies depending on the load characteristics and in 
general can be a non-linear function. 

Utilizing the trapezoidal method to numerically solve (9) 
requires the solution of a set of linear and non-linear equa- 40 
tions: 



The nonlinearity in the above set of equations can be handled 
by the standard Newton-Raphson technique of linearizing 
(11) and iterating until convergence: Equation (10) is a linear 
equation of the form L(x)=0, where we have omitted the 
iteration index k for simplicity. Equation (1 1) is a nonlinear 55 
equation of the form g(x)=0. Let g(x)«G(x) be a linear 
approximation of g(x), linearized around some x=x 0 . Then, 
simultaneously solving L(x)=0 and G(x)=0 yields numerical 
values for x and hence v„. These values are then used to obtain 
a new linear approximation g(x)«G wew (x), and the process is 60 
repeated until convergence. A good choice of the point x 0 for 
the initial linearization at the kth time- step is given by the 
value of v w from the previous time- step. 

A direct implementation of this algorithm requires 
0(pqn 1 3 ) operations , where p is the number of time step s , q i s 65 
the maximum number of Newton-Raphson iterations in each 
time step, and n x =3N. 


We begin by decomposing C, A, and A z as: 



Here C vv denotes the sub-matrix of the capacitance matrix 
that changes amid the simulation, while all other sub -matrices 
remain constant. The matrix C vv captures the drain, gate and 
bulk capacitances of all devices, which are voltage-depen- 
dent, while C cc , and are the capacitance matrices that 
arise from interconnects and are hence constant. 

For typical interconnect structures, the above decomposi- 
tion allows us to manipulate the MNA equations (10) and 
( 11 ): 


jl + l + 4 AlPccA ^' ■U~I~4 Al / V'[)iJ+A 1 v* c + 


( 12 ) 


^4lP«Afi(£ +1 +£)-A,P«C cv (vt +1 -v*) + + v*), 


(13) 


; - ^P C cAl($ +i + $)+ jPccAhil^ 1 + /*) - P cc C cv (v^ +1 - vj), 


c vv vj +1 = C vv vJ - + $ ) + 


(14) 


^Af 2 (/f 1 + Jj) - C vc (v* +1 - vj) + j(I k v +l + I k v ), 


/* +1 = /(v* v +1 ). (15) 

Here r denotes the size of interconnect structure connected 
directly to non linear circuit, and given l=N-r we note that 

C cc eR 1 * 2 , C vv eR rxr . 

P cc =C cc _1 is the inverse capacitance matrix, and A is the 
adjacency matrix of the circuit. A is obtained by first adding 
A^ and A I and then removing zero columns (these correspond 
to intermediate nodes, representing the connection of a resis- 
tance to an inductance). 

Thus far, the analysis is similar to that of the linear ele- 
ments structures described above, with the major difference 
being the addition of (14) and (15), which account for the 
nonlinear elements. We will show here how the linear tech- 
niques can be extended to handle the case when nonlinear 
elements are present. 

For future reference, we will call the technique of directly 
solving (12), (13), (14), and (15) as the “Exact-RLP” algo- 
rithm. It can be shown that the computational complexity of 
the Exact-RLP algorithm is 0(l 3 +pq(l 2 +r 3 )). For large VLSI 
interconnect structures we have l»r, reducing the complex- 
ity to 0(l 3 +pq(l 2 )). 

We now turn to the fast solution of equations (12) through 
(15). Recall that the nonlinear equation (1 5) is handled via the 
Newton-Raphson technique. This requires, at each time step, 
linearizing (15) and substituting it into (14). The resulting set 
of linear equations have very specific structure: 

Equations (12) and (13) are of the form Ax=b where A is 
fixed (does not change with the time-step). Moreover, 
A -1 is typically approximately sparse. 

Equation (14) (after the substitution of the linearized (15)) 
is again of the form Ax=b, where the matrix A is 
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obtained by adding C vv and the coefficient of the first- 
order terms in the linearized equation (15). Recall that 
the matrix C vv captures the drain, gate and bulk capaci- 
tances of all devices. It also contains the interconnect 
coupling capacitances between gates and drains of dif- 5 
ferent non-linear devices in the circuit. As each non- 
linear device is connected to only a few nodes and the 
capacitive effects of interconnects are localized, the A 
matrix is observed to be sparse in practice. Note that A 
changes with each Newton-Raphson iteration and with 10 
the time-step. 

Thus the key computational problem is the solution of a 
sparse time-varying set of linear equations, coupled with a 
large fixed system of linear equations Ax=b with A -1 being 15 
sparse. 

Krylov subspace methods have been shown to work 
extremely well for sparse time- varying linear equations. Spe- 
cifically, the GMRES (Generalized Minimum Residual) 
method of Saad and Schultz allows the efficient solution of a 20 
sparse, possibly non-symmetric, linear system to within a 
pre-specified tolerance. This method performs a directional 
search along the orthogonal Amoldi vectors which span the 
Krylov subspace of A. That is, given an initial guess x 0 and 
corresponding residual r 0 =b-Ax 0 , orthogonal vectors {q 1? 25 
q 2 , . . . , q m } are generated with the property that they span S m , 
the solution search space at iteration m. 


interconnect structure directly connected to non-linear 
devices), by exploiting the Hadamard product formula- 
tion. 

It can also be shown that this formulation facilitates the fast 
matrix- vector products, again in linear time (0(r)), 
which arise while solving linear systems of equations 
with the preconditioner matrix. 

A simple example which best illustrates these advantages 
is a symmetric tridiagonal matrix. 


r ai -bi 
-bi a 2 -&2 


-b n - 2 On - 1 -bn - 1 

i a n , 


(17) 


The inverse of B can be represented compactly as a Hadamard 
product of two matrices, which are defined as follows: 



'«i 

Mi . 

. Mi ' 


f Vi 

Vl ■ 


B~ l = 

Ml 

m 2 . 

. m 2 

o 

V2 

V2 ■ 

■ v n 



m 2 ., 

u n / 


^ v n 

V n ., 



u v 


S m = * 0 + span{ro, Aro, ... , A m r$} (16) 30 

= *o + k(A, ro , m) 

C span{^i, q 2 ... , q m }. 

35 

These vectors are chosen according to the Amoldi itera- 
tion: AQ m =Q„ +1 H„ where Q m ={q x , q 2 ,...,q„} is orthogonal 
and H m eR'" +lx '" is an upper Heisenberg matrix. 

For these methods the choice of a preconditioner matrix M, 
which is an approximation of A, can greatly affect the con- 40 
vergence. A good preconditioner should have the following 
two properties: 

M _1 A«I. 

It must accommodate a fast solution to an equation of the 
form Mz=c for a general c. 45 

FIG. 9 depicts the sparsity structure of the A matrix for a 
circuit example of parallel wires driving a bank of inverters. 

For such a sparsity structure, an appropriate choice of the 
preconditioner could be of the form as shown in FIG. 10. 
Although we have chosen a circuit with only inverters for 50 
simplicity, a more complicated circuit structure would simply 
distribute the entries around the diagonal and off-diagonal 
bands and lead to possibly more off diagonal bands. To see 
this, consider an extreme case where the circuit under con- 
sideration has only non-linear devices and does not comprise 55 
of interconnects. In this case the sparsity pattern of the A 
matrix is as shown in FIG. 11. Therefore, the chosen precon- 
ditioner would encompass not only the sparsity structure 
shown in FIG. 9 but also other sparsity patterns that might 
arise with the analysis of more complicated non-linear 60 
devices. Correspondingly the structure of the preconditioner 
(see FIG. 10) would have additional bands. 

Matrices of the form shown in FIG. 10 have the following 
two properties which make them an ideal choice for precon- 
ditioner. 65 

The inverses of the preconditioner matrix can be computed 
efficiently in linear time, 0(r) (r denotes the size of 


There exists an explicit formula to compute the sequences 
{u},{v} efficiently in 0(n) operations. In this case, if we are 
interested in solving a linear system of equations By=c, we 
only need to concern ourselves with the matrix-vector prod- 
uct B _1 c=y. This computation can also be performed effi- 
ciently in 0(n) computations as outlined below: 


' » (19) 

P «; = Zj u J c j' Pv i = Zj VjCj ' /=1 ’ ■■■ ’ " 
j= i j= i 

yi = Ml P Vl , 

y i = Vi p + m P v . , i = 2, ... , n. 


The above formulation for a tridiagonal matrix could be eas- 
ily extended to handle the more general case when the pre- 
conditioner matrix is a zero padded block tridiagonal matrix 
(matrix with zero diagonals inserted between the main diago- 
nal and the non-zero super-diagonal and sub -diagonal of 
tridiagonal matrix) as in FIG. 10 . Elementary row and column 
block permutations could be performed on such a matrix to 
reduce it into a block tridiagonal matrix. This has been shown 
with a small example as below. 


Ml 

0 

-bi 

0 ' 


(20) and (21) 

0 

«2 

0 

—b 2 



~b\ 

0 

a 3 

0 



0 

-b 2 

0 

04 , 



' Ml 

-bi 

0 

0 ' 



-h 

a 2 

0 

0 

P T , 


0 

0 

a 3 

-b 2 



, o 

0 

—b 2 

a 4 / 




X 
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-continued 

where 

(10 0 0 

0 0 10 
P = 

0 10 0 
,0 0 0 1 

Hence, 



0 Mi 0 Vi 0 V3 0 'j (22) 

0 M2 0 M2 0 V2 0 V4 

o 

mi 0 M3 0 V3 0 V3 0 

0 M2 0 M4 0 V4 0 V4 , 

We have not included block matrices for simplicity of pre- 
sentation, however the zero padded block tridiagonal case is 
a natural extension of the above example. All the entries in 
U,V matrices have to be now replaced by blocks and accord- 
ingly the row and column permutations would be replaced by 
their block counterparts with an identity matrix of appropriate 
size replacing the ‘ones’ in the P matrix. Table 4 gives the 
comparison for Incomplete-LU preconditioner and the zero- 
padded (Z-Pad) preconditioner. Simulations were done on 
circuits consisting of busses with parallel conductors driving 
bank of inverters. ‘Size’ denotes the number of non-linear 
devices. All the results are reported as a ratio of run-time and 
iteration-count (number of iterations for the solution to con- 
verge to within a tolerance of le-10) of Z-Pad to the Incom- 
plete-LU preconditioner. As can be seen from Table 4, Z-Pad 
offers a substantial improvement in run time as compared to 
the Incomplete-LU preconditioner. 


TABLE 4 


Preconditioner comparison. 




Size 



400 

800 

1600 

3200 

Runtime 

.44 

.42 

.42 

.43 

Iterations 

5/10 

5/10 

5/10 

5/10 


We now turn to the solution of equations (12) and (13). As 
mentioned earlier, these equations reduce to the form Ax=b 
with a constant, approximately sparse A -1 . A (corresponding 
to X in (1 2) is composed of L, R and P. Each of these matrices 
has a sparse inverse for typical VLSI interconnects which 
then leads to a approximately sparse A -1 (Note that this argu- 
ment is used for motivating the sparsity inherent in A -1 and 
cannot be used as a theoretical proof for the same). In addition 
this sparsity has a regular pattern which can be explained on 
the basis of how inductance and capacitance matrices are 
extracted. The distributed RLC effects of VLSI interconnects 
can be modeled by dividing conductors into small subsets of 
segments, each of which are aligned. Each of these subsets 
leads to a sparsity pattern (corresponding to a band in A -1 ). 
All the effects when summed up lead to a A -1 matrix that has 
a regular sparsity pattern. Window selection algorithm can 
then be employed to find out the sparsity pattern in A -1 . It has 
been recognized in earlier work that this property (sparsity) 
yields enormous computational savings; it has been shown 
that an approximate implementation of the Exact-RLP algo- 
rithm, referred to simply as the “RLP algorithm” provides an 
order-of-magnitude in computational savings with little sac- 
rifice in simulation accuracy. 


Although X is a dense matrix, X -1 turns out to be an approxi- 
mately sparse matrix. Moreover the matrices X -1 Y, X _1 A l5 
X _1 A 1 P C£ A zl r , X _1 A 1 P cc C (rv are also approximately sparse. 
This information can be used to reduce the computation sig- 
nificantly by noting that each step of trapezoidal integration 
now requires only sparse vector multiplications. Solving 
sparse (23) and (24) along with (15) and (1 6) is termed as the 
RLP algorithm (SASIMI). To analyze the computational sav- 
ing of the approximate algorithm over the Exact-RLP algo- 
rithm, we denote “sparsity index” of a matrix A as ratio of the 
number of entries of A with absolute value less than € to the 
total number of entries. The computation required for each 
iteration of (23) and (24) is then 0((l-v)l 2 ), where v is the 
minimum of the sparsity indices the matrices X -1 Y, X -1 A l5 
X _1 A 1 P cc A zl r , X -1 A 1 P cc C t7V . FIG. 12 provides the average 
sparsity for the matrices for a system with parallel conductors 
driving a bank of inverters. The sizes in consideration are 100, 
200, 500 and 1 000. On top of this the computation time of can 
be reduced to 0(1) by using the windowing techniques (details 
in [16]). Hence the computational complexity of RLP is 0(pq 
(l-v)l 2 ) as compared to 0(pqn : 3 ) for the MNA approach. 

We implemented the Exact-RLP and RLP (SASIMI) algo- 
rithms in C++. A commercially available version of SPICE 
with significant speed-up over the public -domain SPICE has 
been used for reporting all results with SPICE. Simulations 
were done on circuits consisting of busses with parallel con- 
ductors driving bank of inverters, with wires of length 1 nun, 
cross section 1 pmxl pm, and with a wire separation of 1 pm. 
A periodic IV square wave with rise and fall times of 6 ps 
each was applied to the first signal with a time period of 240 
ps. All the other lines were assumed to be quiet. For each wire, 
the drive resistance was 10Q. A time step of 0. 1 5 ps was taken 
and the simulation was performed over 30 ps (or 200 time 
steps). For the inverters the W/L ratio of NMOS and PMOS 
were taken to be 0.42 pm/0.25 pm and 1.26 pm/0.25 pm 
respectively. 

In order to explore the effect of the number of non-linear 
elements relative to the total, three cases were considered. 
With p denoting the ratio of the number of linear elements to 
that of non-linear elements, the experiments were performed 
for p equaling 5, 20 and 50. The number of linear elements in 
the following results is denoted by a. 

We first present results comparing the accuracy in simulat- 
ing the voltage waveforms at the far end of the first line (after 
the inverter load). The RMSE is again use for comparing the 
simulations, defined here as 

S’* 

where v and v denote the waveforms obtained from Exact- 
RLP and SASIMI respectively. 


15 


20 


25 


30 


35 


40 


45 


50 


55 


60 


65 


B [ = PXT [ P t = 
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Table 5 presents a summary of the results from the study of 
simulation accuracy. It can be seen that the simulation accu- 
racy of the Exact-RLP algorithm is almost identical to that of 
SPICE, while the SASIMI has a marginally inferior perfor- 
mance as measured by the RMSE. The error values for 5 
SASIMI are compared simply with the Exact-RLP as it had 
the same accuracy as SPICE results for all the experiments 
run. A plot of the voltage waveforms at the far end of the 
active line, obtained from SPICE, Exact-RLP and SASIMI 
algorithms, is shown in FIG. 13. (The number of conductors to 
in this simulation example is 200.) There is almost no detect- 
able simulation error between the SASIMI, Exact-RLP and 
SPICE waveforms over 200 time steps. To give a better pic- 
ture, the accuracy results reported are for a larger simulation 
time of 2200 time steps. 15 


TABLE 5 


RMSE comparison. 

CT 

p = 5 

p = 20 

o 

Cl 

100 

.0054 

.0053 

.0088 

200 

.0078 

.0052 

.0071 

500 

.0006 

.0022 

.0001 

1000 

.0003 

.0005 

.0004 

2000 

.0003 

.0004 

.0004 


The Hamiltonian of a valley b for electrons, associated with 
the device under consideration, is as follows: 


H b {r) = 


ft d ( 1 d \ d ( 1 d d ( 1 d \ 

2 dx\m b dx) dy\ m b y dy dz\m b dz ) 


+ V(r), 


(25) 


where (m x 6 , m^, m z b ) are the components of effective mass in 
valley b. The equation of motion for the retarded Green’s 
function (G r ) and less-than Green’s function (G^ are: 



-Hbin) \G r b (n,r 2 ,k z ,E)- 


(26) 



(n, r 2 , k z , E)G r b {n, r 2 , k z , E) = <5(n - r 2 ), 


frk 2 

- 2 “Wn) G r b (n,r 2 ,k z ,E)~ 



(n, r 2 , k z , E)G b (n, r 2 , k z , E) = 


(27) 



(r l ,r 2 ,k z ,E)G a b (r l ,r 2 ,k i ,E). 


We now turn to a comparison of the computational require- 
ments between Exact-RLP, SASIMI and SPICE. Table 6 sum- 
marizes the findings. For a fair comparison, our total simula- 
tion time i s compared against the transient simulation time for 3 0 
SPICE (i.e we have not included any of the error check or set 
up time for SPICE). As can be seen from the table, SASIMI 
outperforms the Exact-RLP algorithm and SPICE. For the 
case of 500 conductors with p=50, the Exact-RLP algorithm 
is 390 times as fast compared to SPICE. SASIMI is about 35 
1400 times faster as compared to SPICE, and more than three 
times faster than Exact-RLP. As can be seen, the computa- 
tional savings increase as the ratio of linear to non-linear 
elements is increased from 5 to 50. The savings also increase 
with increase in the size of the problem considered. The 40 
computational efficiency of the SASIMI can be explained on 
the use of sparsity-aware algorithms for both the linear and 
non-linear parts of the problem. 


Given G r and G K , the density of states and the charge density 
can be written as a sum of contributions from the individual 
valleys, 

N(r, k z , E) = V N b (r, k z , E) = -ilm[GJ(r, r, k z , £)], 

b ^ 

p(r, k z ,E) = Yj P^ r ' ('< r - £)]• (29) 


The self-consistent solution of the Green’s function is 
often the most time intensive step in the simulation of electron 
density. It was shown by Svizhenko et al. in “Two-dimen- 
sional quantum mechanical modeling of nanotransistors,” 
Journal of Applied Physics, 91(4):2343-2354, 2002, hereby 


TABLE 6 


Run time (in seconds) comparisons. 



P = 5 



p = 20 



p = 50 


a 

SPICE 

Exact-RLP 

SASIMI 

SPICE 

Exact-RLP 

SASIMI 

SPICE 

Exact-RLP 

SASIMI 

100 

11.96 

1.34 

1.26 

13.73 

.27 

.21 

13.54 

.15 

.12 

200 

100.25 

3.28 

2.68 

68.72 

.64 

.28 

67.68 

.55 

.22 

500 

3590.12 

17.13 

4.872 

1919.21 

13.47 

3.01 

1790.67 

4.58 

1.30 

1000 

>12 hrs 

87.75 

22.71 

>10 hrs 

79.07 

16.49 

>10 hrs 

77.56 

15.20 

2000 

>1 day 

545.6 

78.06 

>1 day 

526.23 

59.33 

>1 day 

408.54 

56.05 


The existing methods for finding the inverse of a block 
tridiagonal matrix suffer from being either numerically 
unstable or heavily memory intensive and hence impractical 
for problems of very large size (e.g. 10 6 xl0 6 ). 

Consider the two-dimensional model of a nano- scale tran- 
sistor, shown in FIG. 14. The body of the transistor is pro- 
jected onto a two-dimensional non-uniform spatial grid of 
dimension N x xN^, where N x (Ny) denote the number of grid 
points along the depth (length) of the device. A brief review of 
the governing physics of this device is provided here. 


incorporated by reference, that the approximate block tridi- 
agonal structure of the left-hand side in (26) and (27) facili- 
tates the efficient calculation of the electron density. Specifi- 
cally, it was demonstrated that only the diagonal entries of G r 
are needed to be computed for such a simulation. 

Svizhenko et al. showed that the problem of computing 
electron densities in a nanotransistor can be reduced to find- 
ing the diagonal blocks of G r , where AG r =I and A is a block 
tridiagonal matrix of the form 


60 


65 
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Ai -Bi 

— B\ A 2 —B 2 

~ B Ny ~ 2 A Ny ~ 1 ~B Ny - 1 

~ B Ny “ 1 A Ny / 


(30) 


where each A z -, B z € C NxX ^ x . Thus A€£ NyNxXNyNx , with N y 
diagonal blocks of size N x each. We will denote A compactly 
as A=trid(A z ,B z ). 

The algorithm in Svizhenko et al. calculates the diagonal 
blocks (D z ) of A -1 in the following manner. 


for large problem sizes. This will preclude it from being 
directly implemented to solve these problems. Alternatively, 
the divide and conquer algorithm described below avoids 
these numerical issues by only handling manageable problem 
sizes at each step. 

In the most simple case, the block tridiagonal matrix, A can 
be decomposed into two sub-matrices connected by what we 
will call a bridge matrix. This concept is illustrated below: 


A = A + XY, 



Gi+rVi+i-BiGfr)- 1 , i= 1, 2,..., N y -l, 

D N =G Ny , 2Q 

A=G,+C? ( 5A.i-B(G < , i=N y -\, N y -2, .... 1. (31) 

The time complexity of this algorithm was shown to be 
©(N^Ny), with a memory requirement of 0(N x 3 N^). 

Alternatively, the inverse of a block tridiagonal matrix can 25 
be computed explicitly by exploiting the Hadamard product 
formulation. When {B z } are non- singular, A is said to be 
proper. In this case there exists two (non-unique) sequences 
of matrices {U Z },{V Z } such that for j>i, (A^X^U^V^. 

Hence, A -1 can be written as 


r Ai -B { 

~b\ A 2 —B2 



~B: 


r A,+i —Bi + 1 

~ B J + 1 A i+ 2 ~ B i+ 2 


S 2 = 


~ B N y ~ 2 A N y ~ I ~ B N y ~ * 
-Bl y - 1 A Ny 


' UiV[ 

ViVl 

UiVl . 

■ v$ y > 


r 0 

0 > 

V 2 V\ ! 

ViVl 

U 2 vl . 

■ VlVSy 

35 

-Bi 

0 

v 3 uf 

V3UI 

ViVl ., 

- UiVty 

X = 

0 

K.— 

cq .. 
I 

VNyVl 

VNyUl 

VNyUf - 

- VNyVJl,, 


, 0 

oj 


40 


The U and V sequences can be efficiently computed in 
0(N_^N x 3 ) operations in the following manner: 

Ui=I, U 2 =B i ~ 1 A l , 

U M =B i - l (A i U r B i _ 1 T U i _ l ), i=2, . . . , N y -1, 

V N y T= { A N y U N y -B Ny _ l T U Ny _ 1 ) 1 , 

V Ny -i T =VN y T A N BN y -i l , 

Vi T =(Vi + i r A i+ i-V i+2 T B i+i T ), i=N y -2 , . . . , 2, 1. (32) 

I denotes identity matrix of appropriate size. 

The matrices encountered in the simulation of nanotrans- 
istors enjoy the property that the diagonal blocks {A z } are 
tridiagonal while the off-diagonal blocks {B z } are diagonal 
(this is due to the tight-binding approximation of the Hamil- 
tonian). The complexity of the above parameterization is then 
reduced to 0(N y N x 2 +N x 3 ). It is worthwhile to note here that 
the complexity of algorithm in (31), which was given to be 
0(N y N x 3 ), does not change based upon these properties (al- 
though the actual run time is reduced). Therefore, the reduced 
complexity of the Hadamard formulation makes it an ideal 
choice for the solution of this problem. However, it has been 
shown that the above recursions can be numerically unstable 


Notice that the block tridiagonal structure is preserved in 
each sub -problem, and they are joined by the N x xN x bridge 
matrix B z . The structure of and S 2 facilitates the use of 
45 Hadamard based methods (32) for determining the solution of 
each sub -problem. However, the question then becomes how 
to determine the global solution from each sub-problem solu- 
tion and their corresponding bridge. This problem can be 
50 resolved by the use of a fundamental result of linear algebra. 
The matrix inversion lemma describes the effect of a low rank 
correction term on the inverse of a matrix, and can be used in 
the following way: 


55 


60 


65 


A~ l =(A + XY) 1 

= a- 1 -{a- 1 x){i+ya- 1 x)-\ya~ 1 ) 

where, 

i-CiBi) 0 
0 (- c 2 Bj ) 



-S^U, 1 )Bf 

/ 


( 33 ) 
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-continued 

o (ciy | 

(Cl) 0 J' 


Here, C 1 =S 1 1 (:,i) and C 2 =S 2 1 (:,1) denotes the last (first) 
block columns of S 1 _1 (S 2 _1 ) respectively. 

Although (33) can be used directly to solve the problem 
described above, its lack of computational efficiency pre- 
cludes its use in this case. However, it is important to note that 
the adjustment term, (I+YA -1 X) -1 , depends only on the cor- 
ner blocks of the inverses of the sub-problems. This observa- 
tion provides the basis for a formulation of the solution for the 
more general case. Specifically, the divide and conquer algo- 
rithm of the present invention addresses the issue of how to 
combine multiple sub-problem solutions in both a memory 
and computationally efficient manner. 

An overview of the divide and conquer algorithm is pro- 
vided here. 

The procedure begins with separating the block tridiagonal 
matrix A into D subproblems, each joined by a bridge matrix. 
The procedure presented previously motivates the formula- 
tion of the global solution by combining the sub-problems in 
a simple radix 2 fashion. However, this method offers no 
improvement in terms of memory consumption and is com- 
putationally intensive. Alternatively, matrix maps are created 
to capture the effect of each combining step without perform- 
ing all of the associated computation. Adjustments to the 
matrix maps at each combining stage are not constant and 
must be modified to parallel the procedure above. These maps 
can then be used in the final stage to transform the subproblem 
solutions into the global solution. 

Each combining step in the divide and conquer algorithm 
will have an associated entry in the “job” list pointing to the 
first and last sub -problems along with the corresponding 
bridge point. For example, (1~2,3~4) describes the action of 
combining joined sub -problems 1 and 2 (S^) to joined sub- 
problems 3 and 4 (S 3 _ 4 ) by use of the bridge matrix between 
problems 2 and 3. The corresponding entry in the job list 
would be of the form: [start, stop, bridge]=[l, 4, 2]. 

This formulation lends itself directly to parallel implemen- 
tation, since computations associated with non-overlapping 
combinations can be farmed out across multiple systems. 

To model the effect of any combining stage in the job list, 
it is necessary to know the comer block elements from the 
inverse of each combined sub -problem. This process can be 
easily illustrated by a simple example. Suppose once again 
that the combination stage is (1~2,3~4), let 
and Q 2 =S 3 _ 4 . The comer block elements of Q : _1 and Q 2 _1 can 
be found using the parameterization given in (32). 
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where n=size (SJ+size (S 2 ), and m=size (S 3 )+size (S 4 ). 

It would be impractical to recalculate the inverse of each 
joined problem for each combining stage. Alternatively, 
1 5 matrix maps are used to efficiently produce the required block 
entries. 

Matrix maps are created to produce the cumulative effect of 
each combining step associated with a particular sub-prob- 
20 lem. There are a total of eight N x xN x matrix maps {M.} for 
each sub-problem. The process of updating matrix maps can 
be broken down into two categories: Adjustments to Upper 
sub -problem and those to Lower sub-problems, the distinc- 
tion being their location with respect to the bridge point. This 
25 procedure will be illustrated using the above example. The 
“adjustment” matrix, ADJ, for the combining step is defined 
as follows: 


30 


Z, = 

Z 2 = -B ml U n V T n , 
P=(/-Z 2 Z,)-‘, 


(34) 


(I + Z l PZ 2 

-ZiP\ 

( ADJn 

ADJ 12 ) 

< -PZ 2 

p ) 

~{adj 2 i 

ADJ 2 2 J 
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The matrix maps associated with the Upper sub-problems (S : 
and S 2 ) are updated in the following manner. 

45 V n T ADJ l2 B n ^)M 3 -^M, 

M 2 (U i V n T ADJ l2 B n+l )M 4 -^M 2 

(O l f r m ADJ 22 B n+l )M 3 -+M 3 

50 (&iK T ADJ 22 B n+l )M 4 -*M 4 

M 5 -M 3 T (ADJ l2 B n+l )M 3 ^M 5 

M 6 -M 3 t (ADJ \ 2 B n+l )M 4 -*M 6 
55 

M 7 -M 4 T (ADJ l2 B„ +l )M 3 ^M 7 


M s -M 4 T (ADJ l2 B n+l )M 4 -^M s (35) 

Those associated with the Lower sub -problems (S 3 and S 4 ) 
60 are updated in the following manner: 

{U { V n T ADJ { 1 B n+1 r )M 1 ^M 1 
(Ui V„ T ADJ n B n+l T )M 2 ^M 2 


M 3+ (V m U 1 T ADJ 2l B n+l T )M l ^M 3 


M 4 +{V m U l T ADJ 2l B n+l T )M 2 ^M 4 
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Ms-Mi T (ADJ 2l B n+ 1 
M 6 -Mi T (ADJ 2l B n+l T )M 2 —>M 6 
M 7 -M 2 T (ADJ 2l B n+l t )M^M 7 

M 8 -M 2 T (ADJ 2l B n+l r )M 2 ^M 8 (36) 

The above procedure for modifying the matrix maps is 
repeated for each entry in the job list. In the final stage the 
matrix maps are used to generate the diagonal blocks of A -1 . 
It is important to note that this scheme fits nicely into a 
parallel framework due to the fact that systems handling 
Upper or Lower sub -problems would only have to trade the 
limited amount of information in (34) to modify the matrix 
maps they are governing. 


Pseudo -Code 

1. For each sub-problem in {1,2,...,D} : 

• Determine corner blocks from inverse of sub -problem 

• Associate bridge matrix (excluding problem D) 

• Initialize matrix maps 

2. Generate list of sub-problem combinations radix 2 

3. Adjust mappings for each combining step 

4. Compute diagonal elements for each division 

5. Apply matrix maps to transform sub-problem solutions to the global 
solution 


The time complexity of the divide and conquer algorithm is 
0(N x “N^+N x 3 D log 2 D) where a is defined to be the order 
associated with a N x xN x matrix multiplication (typically 
a«2.7). The memory consumption is 


The divide and conquer algorithm, along with the algo- 
rithm in Svizhenko et al. have been implemented, in Matlab, 
on a single 32 -bit x86 Linux workstation. All results reported 
are for test cases on a MIT well-tempered 25 nm device-like 
structure. 

Table 7 shows the run-time comparison of the two algo- 
rithms across the cases: N x =100, ^=3,000-80,000. Notice 
that the algorithm given in Svizhenko et al. was unable to 
perform past the case N v =l 1,000 due to memory restrictions. 
However, the divide and conquer algorithm was able to 
handle these cases without encountering memory overflow. 
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FIG. 15 shows the ratio of memory consumption of the 
algorithm in Svizhenko et al. as compared to the divide and 
conquer algorithm for varying number of blocks per division 
(sub-problem size). FIG. 16 shows the same ratio for varying 
5 number of divisions. 

While the invention has been illustrated and described in 
detail in the drawings and foregoing description, the same is 
to be considered as illustrative and not restrictive in character, 
10 it being understood that only the preferred embodiment has 
been shown and described and that all changes and modifi- 
cations that come within the spirit of the invention are desired 
to be protected. 

What is claimed is: 

1 . A system for simulating operation of a VLSI intercon- 
nect structure having capacitive and inductive coupling 
between nodes thereof, comprising: 

20 a processor; and 

a memory, said processor configured to: 

obtaining a matrix X and a matrix Y containing different 
combinations of passive circuit element values for 
25 said interconnect structure, said element values for 

each matrix including inductance L and inverse 
capacitance P, 

obtaining an adjacency matrix A associated with said 
30 interconnect structure, 

storing said matrices X, Y, and A in said memory, and 

performing numerical integration to solve first and sec- 
ond equations each including as a factor product of 
35 inverse said matrix X (X -1 ) and at least one other 

matrix, said first equation including X -1 Y, X -1 A, and 
X -1 P, and said second equation including X -1 A and 
X _1 P. 

2. The system of claim 1, wherein said matrices X and Y 
each contain resistance values in addition to L and P. 

3. The system of claim 2, wherein said first equation is 
substantially of the form, 

45 


jf +1 = X" 1 Y$ + X~ l A^ n + ^ X~ l AP(I k s +l + /*), 
50 


where v„ and i z are node voltages and inductor currents, 
respectively, A is an adjacency matrix for the circuit, and 


TABLE 7 


Run time (min) comparison 






Size = N x * N v (xlO 6 ) 






0.4 

0.5 

0.6 

0.7 

0.8 

0.9 

1.0 

1.1 

2.0 

4.0 

8.0 

Divide and Conquer Algorithm 

5.38 

6.66 

8.10 

9.41 

11.13 

12.25 

13.75 

14.99 

27.59 

58.86 

107.0 

Algorithm in Svizhenko et al. 

1.22 

1.51 

1.82 

2.12 

2.43 

2.72 

3.01 

3.33 

— 

— 

— 
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l s is a current source vector, and wherein said second 
equation is substantially of the form 


X- Y Avi +i : 


X^Av^-X^APA 7 ^ 
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